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We present a comparison of methods for treating the electrostatic interactions of 
finite, isolated systems within periodic boundary conditions (PBCs), within Den- 
sity Functional Theory (DFT), with particular emphasis on linear-scaling (LS) DPT. 
Often, PBCs are not physically realistic but are an unavoidable consequence of the 
choice of basis set and the efficacy of using Pourier transforms to compute the Hartree 
potential. In such cases the effects of PBCs on the calculations need to be avoided, 
so that the results obtained represent the open rather than the periodic boundary. 
The very large systems encountered in LS-DPT make the demands of the super- 
cell approximation for isolated systems more difficult to manage, and we show cases 
where the open boundary (infinite cell) result cannot be obtained from extrapola- 
tion of calculations from periodic cells of increasing size. We discuss, implement and 
test three very different approaches for overcoming or circumventing the effects of 
PBCs: truncation of the Coulomb interaction combined with padding of the simula- 
tion cell, approaches based on the minimum image convention, and the explicit use 
of Open Boundary Conditions (OBCs). We have implemented these approaches in 
the ONETEP LS-DPT program and applied them to a range of systems, including 
a polar nanorod and a protein. We compare their accuracy, complexity, and rate of 
convergence with simulation cell size. We demonstrate that corrective approaches 
within PBCs can acheive the OBC result more efficiently and accurately than pure 
OBC approaches. 
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I. INTRODUCTION 



Density Functional Theory (DFT)P21 is widely and routinely used for computational elec- 
tronic structure simulations due to its favorable balance of speed and accuracy. However, 
making DFT simulations scale well to the numbers of atoms required to study large complex 
systems such as proteins and nanostructures presents significant challenges. Various linear- 
scaling approaches to DFT have emerged over the last two decades to meet this challengd^^. 
Several of these methods use basis sets which are related to plane waves and require periodic 
boundary conditions (PBCs). The plane-wave pseudopotential approach has been developed 
with crystalline systems in mind, and as these are genuinely periodic, the treatment of elec- 
trostatics in the framework of PBCs was a natural choice with significant advantages. In 
reciprocal space, the Hartree interaction is diagonal, so the Hartree potential and energy are 
easily obtained using Fast Fourier Transforms (FFTs). Furthermore, the plane- wave basis 
set is systematic in the sense that it provides a uniform description of space and can be 
improved by increasing the value of one parameter. 

However, the increasing use of linear-scaling DFT (LS-DFT) in large systems highlights 
long-standing issues in electronic structure methods relating to the treatment of electrostatic 
interactions, i.e. the long-ranged parts of the Coulomb interaction between electron density 
and electron density ('Hartree' terms), electron density and ion cores, and between ion cores, 
under PBCs. 

Bulk systems can be genuinely periodic and then the infiuence of periodic replicas is 
desired; however, to allow simulation of finite, isolated systems within PBCs, the supercell 
approximation is widely usecP^'^. This involves the replacement of a genuinely isolated 
system with a lattice of periodic replicas, with vacuum 'padding' surrounding the system to 
reduce the infiuence of periodic replicas on each other. While this is a reasonable approach, 
it introduces finite size errors whereby the total energy varies with supercell size. 

The use of a supercell is frequently a well-controlled approximation: that is to say, by 
increasing the size of the cell and thus the distance between periodic images, one rapidly 
approaches the true isolated, non-periodic limit. For example, in the case of relatively small, 
charge- neutral molecules without significant dipole moment, one needs to ensure simply that 
the charge densities of periodic replicas do not overlap to any significant extent. In other 
cases, the amount of vacuum padding required to reach this limit can become prohibitively 
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large. The slow decay of the interaction of periodic replicas of a monopole charge, as l/i?, 
means that the infinite limit is impossible to reach in practice for charged systems. Similarly, 
for highly-elongated charge-neutral systems possessing a large dipole moment (such as in a 
simulation of a polar semiconductor nanorod), the simulation cell would likewise need to be 
unfeasibly large to prevent interactions between periodic images of adjacent rods. Clearly 
the isolated limit cannot always be found simply by extrapolating to infinite supercell size. 
This issue is exacerbated as the isolated molecules and their dipole moments become larger. 

To address this problem, a large range of techniques that aim to either reduce or eliminate 
the effects of the PBCs on the electrostatics of grid-based electronic structure calculations 
have been developed over the recent year^^^'^. These include methods which attempt to 
formulate an a posteriori correction term to add to the energ}^* ^^ * ^^ ! on the basis of a multi- 
pole expansion of the localised charge, having first inserted a uniform periodic background to 
counter any monopole charg^^, methods which formulate a more complex form of 'counter- 
charge' which counteracts the periodic interaction j^^ * ^^ * ^^ * ^^ * ^° l, and methods that modify the 
form of the interaction in real or reciprocal space in order to avoid the existence of periodic 
interactions in the first plac^ZEQEam. 

In this paper, we examine, implement and compare three different approaches fulfilling 
these criteria: truncation of the Coulomb interaction in real space, referred to here as 'Cutoff 
Coulomb' (CC)llI12il the approaches of Martyna and Tuckerman (MT) and Genovese et al, 
which replace the periodic Coulomb interaction with a Minimum Image Convention (MIC) 
approach to the Coulomb potentiaP^, and the replacement of PBCs with Open Boundary 
Conditions (OBCs) using a multigrid approach to the Poisson equationpHill These methods 
are implemented and tested on a range of systems representing typical cases with challenging 
electrostatic properties. We compare their accuracy, convergence properties, complexity and 
computational overhead, and summarise the advantages and disadvantages of each. 

Throughout this work, we employ linear-scaling DFT with the ONETEP codd^Sl^ 
while our findings will be applicable to all electronic structure methods, linear-scaling or 
otherwise, we focus in particular on the challenges encountered applying these methods 
to large, complex systems. System size can be measured either by the number of atoms 
N included in the simulation, or by the volume V of the simulation cell — the latter 
being particularly relevant in the case of isolated systems. ONETEP combines linear-scaling 
computational effort, in that the total computational time for a simulation of N atoms can 
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be made to scale as 0{N), with near-independence of the computational effort on the amount 
of vacuum padding (i.e. nearly independent of V at fixed A^), and systematic control of the 
accuracy with respect to the basis, akin to that of plane-wave DFT. The requirements on any 
method used to treat electrostatic interactions are therefore that it must have systematically 
controllable accuracy, must not impose too high a computational overhead, and must have 
low-order scaling with and V. 

II. ELECTROSTATICS IN LINEAR-SCALING DENSITY FUNCTIONAL 
THEORY 

The calculations in this work are performed with the ONETEP Linear-Scaling DFT 
approach. Like most linear-scaling approaches to DFT, ONETEP uses the density matrix 
rather than eigenstates of the Hamiltonian, representing the single-electron density matrix 
p(r,r') in terms of nonorthogonal localised orbitals 0Q,(r) and a 'density kernel' K'^^ as 

p(r,r') = 0.(r)ir"%(r'). (1) 

The Einstein convention of summation over repeated Greek indices will be employed through- 
out. Using the density matrix, the electron density n(r) can be found from 

n(r) =p(r,r) = 0„(r)ir°%(r) . (2) 

Where ONETEP differs from most linear-scaling approaches is that the local orbitals, re- 
ferred to as Nonorthogonal Generalised Wannier Functions (NGWFs)'2S, are themselves ex- 
pressed in a systematic underyling basis of periodic-sinc functions (psincs), and are therefore 
systematically convergeable. This is achieved by a double-loop optimisationp3 of both the 
coefficients Cia of the psinc functions -Di(r) describing each NGWF and the elements of the 
density kernel i^'"^: 

= mm LUQ.}) , (3) 

where L represents optimisation with respect to the density kernel, a generalisation of the 
occupancies, through: 

L{{Q^})= mmE{{K^^};{a^}). (4) 

{_ft-Q/3} 

This results in a method with controllable accuracy and systematic convergence of total 
energies and forces with respect to basis size, equivalent to the plane-wave approachP'^, 
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in systems of tens of thousands of atomd^214i|_ Convergence is controlled by varying the 



spacing of the psinc grid, in a manner equivalent to varying a plane-wave cutoff, described 
by a cutoff energy E^nt, and by varying the cutoff radii of the spherically-truncated NGWFs, 
described by a sphere radius R^. To achieve true asymptotically linear scaling behaviour, 
it is also necessary to truncate the range of the density kernel K"^ so that elements for 
NGWFs centred on distant atoms for which |Ra — R/s] > -Rk are set to zero. However, this 
latter form of truncation is only necessary in very large systems and will not be considered 
in this work. 

This accurate and systematic approach to linear-scaling total energy calculations demands 
that all aspects of the calculation be carried out with high accuracy, including the long- 
range electrostatic part. The electrostatic energy comprises the Hartree term, -ErM, which 
is the classical density-density interaction; the local pseudopotential term, E\ocps[n], which 
is the interaction of the electron density with the long-ranged part of the potential resulting 
from the ion cores; and the interaction between the ion cores, -Eion-ion- It should be noted 
that during the optimisation of the kernel and NGWF coefficients i^"^ and Cia, the full 
interacting energy is minimised by conjugate gradients process, meaning that no mixing of 
densities is required at any point. The problem, then, becomes one simply of evaluating 
i^nM and Vh['^](i") for a given density n{r) (which always integrates fo the number of 
electrons A^^e). 

To be absolutely clear on the formalism involved, we will briefly re-visit the standard 
approach, making careful distinctions on how the expressions and their meaning vary under 
PBCs and under OBCs, where the potentials tend to zero at infinity. In both cases, the 
Hartree energy can be obtained as E-^ = | / n{r)Vu{r) dr, where the Hartree potential VH(r) 
resulting from a density n{r), is formally obtained by solving the Poisson equation: 



Note that we are working in atomic units, for which I/eq = An. This can in general be solved 
through the use of the corresponding Green function G(r, r') = — l/47r|r — r'|, producing 



This result builds in the OBC definition that the potential goes to zero at infinity, and 
cannot be used directly to evaluate E^ or VH(r) under PBCs as the integral has infinite 
value at all r for periodic n(r'). 



VVH(r) = -Ann{r) . 



(5) 
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When PBCs are used Eq. ([s]) is only valid for charge distributions of zero charge per 
simulation cell. If the total charge on one cell Q = Jq n{r) dr is non-zero, Eq. (j5| is modified 
to the following form: 

VVH(r) = -47r(n(r) - q/Q) , (6) 

where Q is the volume of the simulation cell. This is equivalent to the insertion of a uniform 
background charge density of equal and opposite charge to n{r) so that the total charge is 
zero. A periodic density will result in a periodic potential and in this case we can re-write 
both sides of Eq. ^ in terms of their discrete Fourier transforms and rearrange to obtain 

VniG) = — {h{G) - q5a,o) . (7) 

Note that Eq. ([T]) makes clear the utility of a reciprocal space approach to calculating 
Vh(G), even outside of a genuinely periodic situation: the Coulomb interaction is diagonal 
in reciprocal space, so Vh(G) can be obtained trivially from n{G). After obtaining ^^(r) 
by an inverse FFT, the integral E'h = ^ ^(i")^h(i") c^r can be performed only over one 
simulation cell to obtain the Hartree energy per simulation cell. 

In PBCs the potential is, by definition, the result of contributions from not just the n{r) of 
the home simulation cell but also from the densities of an infinite number of periodic replicas 
of that cell. A periodic function that can be constructed in this way is demonstrated with the 
example at the top panel of Figure [T| As we have already mentioned, the potential and the 
electrostatic energy diverge for non-zero total charge in the simulation cell (or equivalently 
when n{G = 0) is nonzero). To avoid this divergence one must set n{G = 0) to zero for 
each component making up total charge density (including the ion charges) to ensure that 
the result is finite. Having made this choice however, one alters the problem being studied 
as the potential Vii{r) obtained is that resulting not just from the infinite periodic array of 
n{r), but also from a neutralising charge distribution, which is usually taken to be a uniform 
background charge over the whole cell. 

The same arguments apply to the other electrostatic terms, by replacing the electron 
density n{r) with the charge density of the ions, in the form of a collection of point charges. 
For an isolated system, the energy of interaction of the ions is of course simply 
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Figure 1 . Different ways of making a function obey periodic boundary conditions inside a simulation 
cell, demonstrated for a Gaussian function. Top panel: The Fourier transform approach. The 
resulting function is the same as the one that would be obtained by a superposition (sum) of 
periodically repeated Gaussians. Bottom panel: The Minimum Image Convention (MIC) approach: 
the resulting function is the same as the one that would be obtained by having a single Gaussian 
in the simulation cell and making it periodic by applying the MIC. 

while under PBCs, in the presence of the neutrahsing background, the energy of interaction 
per unit cell is most commonly calculated using the Ewald techniqu^^. 

The influence of periodic neighbours will affect (polarise) the charge distribution during 
a self-consistent electronic structure calculation. Therefore, it should be immediately clear 
that no a posteriori approach to correcting total energies obtained from a simulation under 
PBCs can be completely successful in providing total energies that match those of an isolated 
system as even after the "removal" of the periodicity the density will remain distorted to 
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what it was in the periodic calculation. Here we examine three approaches that are applied 
within the self-consistent procedure and therefore are able to correct not only the energy 
but also the potential. 



III. CUTOFF COULOMB INTERACTIONS 

One way to avoid the effects of PBCs which are intrinsic to the discrete Fourier representa- 
tion of the Coulomb potential is to use a modified form for the Coulomb potential. One such 
possibility is the use of a "cutoff" form of the Coulomb interaction. This allows the usual 
Fourier transform-based approach to be used, including a nominally periodic simulation cell, 
but truncates the Coulomb potential so that it is confined within the primary simulation 
cell. The approach has been applied by several previous work^^^'^ll a^fj jg implemented in 
several codep3E31. 

The essence of the cutoff Coulomb approach is that the periodic, background- neutralised 
Coulomb potential VEw(r) is replaced with the bare Coulomb interaction, truncated so as to 
prevent any part of the simulation cell feeling the potential from any neighbouring copy. This 
removes the need for the canceling background, even though the charge density is periodically 
repeated. Some new complications arise however as the cutoff Coulomb potential needs to 
be generated in reciprocal space. 

To retain the simplicity of having an interaction that is diagonal in reciprocal space, but 
still avoid the influence of periodic replicas, one can use the following form for the Coulomb 
potential 



TZi is a region of a size and shape chosen such that when centered at any point r at which 
Vii{r) is required (this may be anywhere inside the main simulation cell, or it may just be 
anywhere where the density is nonzero), TZi encloses all r' for which n{r + r') 7^ 0. The 
Hartree potential is now obtained as the convolution of the cut-off Coulomb operator and 
the density 



Jn 

The simplest shape for TZi is a sphere of radius Rc, for which Vq^^^^^t) = 0(|r| — i?c)/|r| 
where is the Heaviside step function. In this case, the Fourier transform of the interaction 






(10) 
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is well-known: 

V^-^-^(G) = 'Jgi^"^^^^ • (11) 
As this function does not have a singularity at G = the Hartree potential is obtained in 
reciprocal space as its product with n{G) as in Eq. ([T]) but without the q term as there is 
no longer the need to include a uniform background charge. A spherical cutoff removes the 
periodicity in all three spatial dimensions. If periodicity is retained in one or two dimensions 
there are corresponding forms for Vcc(G) to account for these wire (ID periodicity) and slab 
(2D periodicity) geometries. A comprehensive study was made by Rozzi et. al^ describing 
the terms of the cutoff Coulomb interaction for each geometry. 

In a practical calculation, the electron density n{r) on a real space grid over the original 
simulation cell is transferred to a grid for a larger 'padded' cell of size Lpad and padded with 
zeros, then Fourier transformed to give npad(G). The terms of Vcc(G) are calculated for this 
reciprocal space grid in advance and stored, and are used to multiply the Fourier components 
^pad(G) whenever the Hartree potential is required. Reverse Fourier transforming these 
components gives VH,pad(r) from which the values of lii(r) on the original cell are extracted. 

The corresponding cut-off form of the Coulomb interaction must also be used in place 
of the long-ranged Coulombic tail of the ion cores in the local pseudopotential V[ocps(i")- To 
achieve this, Viocps(G) is calculated over the whole padded grid, replacing the j^Zion term 
by Vcc(G)Zion for the relevant form of cutoff Coulomb interaction. This is then transformed 
to real space by Fourier transform and extracted to the standard grid to give the required 
form of V|ocps(r). Similarly, the periodic Coulomb and Ewald terms in the calculation of the 
forces acting on the ion cores are replaced by their cutoff Coulomb forms. 

The computational overhead of the method during an SCF calculation compared to the 
traditional PBC Fourier transform Coulomb approach consists of three parts: transfer of the 
calculated density from the original grid to a larger, padded grid, calculation of the forward 
and backwards Fast Fourier Transforms required for the Hartree potential on the larger 
grid, and extraction of the calculated potential from the larger grid back to the original 
one. The first and last of these are in general comparatively trivial and take very little 
time. Performing the FFT on the larger grid often incurs a considerable slowdown relative 
to performing it on the original grid, but nevertheless, generally speaking, this part of the 
calculation takes a almost negligible fraction of the total computational time for large enough 
systems. 
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Figure 2. Illustration of the cell sizes Lcelb -^pad ^-nd cutoff radius i?c required for the spherical 
cutoff Coulomb approach. Rc must be at least as large as the largest distance between any two 
non-zero charges in the system (this is trivially satisfied if Rc > \/3iccll)- In order for the periodic 
densities not to impinge on each other, Lpad ^ (-^mol + ^c) rnust be satisfied, where Lmol is the 
extent of the system (again, defined as maximum distance between two non-zero charges) in any 
Cartesian direction. 

When simulating an isolated object such as a nanocrystal or nanotube with a high aspect 
ratio, the geometry of the system requires that we use a simulation cell that is very long 
in one dimension (the x direction here) and comparatively small in the other two {y and 
z). Performing cutoff Coulomb calculations with a spherical cutoff would rapidly become 
impractical as the length of the system is increased, since for a sphere geometry, we would be 
required to embed the original cell in a padded cell with all the side lengths , Ly, > Re- 
in such cases, we need to define a geometry for the cutoff Coulomb interaction such that 
the cutoff range can be very long in one direction and shorter in the other two. One obvious 
choice for a long, thin system is to cut off the Coulomb interaction on the surface of a 
cylinder. In this case, the integrals required to evaluate the coefficients are not analytically 
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solvable but can be put in a form amenable to numerical evaluation. Appendix A gives 
details on the evaluation of the Fourier coefficients of the interaction for a cylindrical cutoff. 
With an efficient system for evaluating the terms Vcc(G) numerically, the interaction can 
be calculated rapidly in advance and reused, and simulations of isolated high aspect ratio 
systems can proceed within cells of feasible size. 



IV. MINIMUM IMAGE CONVENTION 

An alternative technique for avoiding periodic interactions is the class of approaches which 
includes those of Martyna-TuckermarP^l and Genovese et aP^^. The essence of these, which 
we will call Minimum Image Convention (MIC) approaches is that the form of the Coulomb 
operator is modified in a way that is still periodic (as this is unavoidable if standard FFTs 
are to be used) but which nevertheless removes contributions from neighbouring cells. 

To see how this is achieved, we consider first the Fourier transform of a function /(r), 
defined as 

/(G) = [ e-^^- /(r) dr. (12) 

Jail space 

In PBCs, a discrete set of wave vectors G are used to expand functions in Fourier space. 
These wavevectors are chosen by the requirement that they need to be commensurate with 
the simulation cell. Therefore, given the expression for /(G), the real space representation 
of the function /(r) under PBCs is the following: 

/per(r)= J2 /(G)e^^-^ . (13) 
G e cell 

This is an exact result and shows that the Fourier representation of /(r) in the simulation 
cell is a periodic function /pcr(i") with the periodicity of the simulation cell. It is important to 
notice is that this function is constructed as a superposition of periodically repeated functions 
/(r), one in each cell. This is demonstrated for the example of a Gaussian function in the top 
panel of Figure [T| where its resulting periodic form in one simulation cell is provided, as it 



would be generated in real space as a Fourier expansion by Eq. (13). This result implies that 
periodic interactions are unavoidable if the potential is constructed by approaches based on 
Fourier transforms in the standard simulation cell, as PBCs are implicit in such procedures. 
However, MIC approaches are designed to avoid the part of the Coulomb interaction which 
produces this undesired long-ranged interaction. 
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We have implemented the Martyna-Tuckerman approachP^, in which the Fourier method 
is used to construct not the periodic function /pcr(i') but the periodic function /Mic(r) 
which results by making /(r) periodic over a single simulation cell using the MIG^^. A 
similar approach can also be employed in Quantum Monte Carlo calculations, via the "Model 
Periodic Coulomb" approachP^'^. The distinction between fperi^) and /Mic(r) is clarified in 
Figure [l] where the bottom panel demonstrates the construction of /Mic(r) for the example 
of a Gaussian function. 

To work with this formalism we need to determine the Fourier transform /(G) that will 
produce the desired /mic(i") 

/Mic(r) = 7(G)e'''-^ • (14) 

G cell 

As this method is intended for dealing with the Coulomb potential, from now on we will 
fix the function /(r) to be equal to 0(r) = 1/r so that we can focus on particular issues 
that arise in this case. In determining the form of 0(G) we need to deal with the extra 
complication of the singularity of the potential at r — > (short range) and at G — > (long 
range). The Coulomb potential is partitioned as 

1 erffo^T') erfcfo^r) 

- = + — = 01ong(r) + 0short(r) , (15) 

ly lyt lyi 

where a is a convergence parameter which determines the region where the transition from 
short to long-range terms takes place. Assuming that the simulation cell is large enough so 
that 0short(G) — 0short(G) ouly the long range form 0iong(G) needs to be determined. The 
desired Fourier transform is expressed as 

0(G)=01ong(G)+<^short(G) (16) 
= [^long(G) - 01ong(G)] + [01ong(G) + 0short(G)] 

= 0screen(G)+0(G) , (17) 



where the explicit expression for 0short(G) is 

0short(G) = — 



exp 



(18) 



Eq. (16) can also be further expanded to the form shown in Eq. (17) which demon- 



strates that the MT formalism is equivalent to augmenting 0(G) with a "screening potential" 
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(G) which cuts off the interactions from the periodic images of the simulation celL In 



practice, we compute 0(G) according to Eq. (16) and we distinguish two cases: G 7^ and 
G = 0, which must be treated separately. 



The function 0iong(G) for G 7^ is obtained as 



long 



'G) 



n r 



(19) 



where the above integral is computed as a sum over the simulation cell grid points as this is 
an exact expression for the wavevectors G which are commensurate with the simulation cell. 
The above expression is the desired one as the term erf (ar) / r does not contain contributions 
from periodic images. It also does not contain a singularity at r = so the evaluation of 
this integral poses no difficulties. The complete expression for 0(G) is obtained as the sum 



of the terms Eq. (18) and Eq. (19). 



To find the G = term, we need to consider the limit of Eq. (18) as G goes to zero 



''short 



lim ( 

G->0 



lim — 

G^O 
71 



'G) 



1-1 



+ 



G^ 



4a2 8 



+ 



(20) 



and taking this into account, Eq. (16) becomes 



0(0) 



lone 



(0) + 



•'short 



0) 



erf(ar) , vr 



(21) 



where the integral in the above expression is again evaluated as a sum over the simulation 
cell grid points as the integrand does not contain a singularity at r = 0. 

In order to use the MT potential in practical calculations, we need to ensure that appro- 
priate conditions are obeyed as regards the relative sizes of the simulated molecule and the 
simulation cell. From the example in the bottom panel of Figure [T] we can see that the length 
that a simulation cell can have in any direction needs to be at least twice the length of the 
molecule being simulated. In the opposite case unphysical interactions will be introduced as 
some charges on the molecule will be experiencing the Coulomb potential from other parts 
of the molecule (as they should) while other charges will experience the potential from a 
periodic image (which they should not). 
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In our implementation the Hartree potential is generated in reciprocal space from the 
electronic density as a product with the Fourier transform of MT potential 0(G) 

Fh(G) = 0(G)n(G) . (22) 

In a similar way, the local pseudopotential is obtained in reciprocal space as a sum of short 
and long range terms 

^locps(G) = V^locps,short(G) + V^locps,long(G) . (23) 

For an ion with charge — Z, (following the established electronic structure theory convention 
of taking the ionic potential as negative), the periodic Coulomb component is subtracted 
from the pseudopotential to obtain its short range part 

^locps, short 

(G) = yiocps(G) + Zci>{G) (24) 

and the long range part is obtained as the MIC Coulomb interaction 

Flocps,long(G) = -^0io„g(G) . (25) 

Finally the core-core interaction energy is obtained as a Coulombic sum between point charge 
interactions in the simulation cell according to Eq. (|8|. 

Genovese et al^^^ proposed an approach that is rather similar in principle but in prac- 
tice has some different properties. They described a wavelet-based approach to calculating 
the MIC Coulomb interaction. The charge density is expanded using interpolating scaling 
functions^'' of order m (typically m = 14). This guarantees that when a known continuous 
charge density is represented, the first m moments are preserved. Although most practical 
methods do not attempt to represent given continuous charge densities, this approach is 
useful when using pseudopotentials of the form proposed by Goedecker et al^. The repre- 
sentation of the Coulomb operator is made separable by employing an expansion in terms of 
Gaussian^. The resulting one-dimensional integrals can be calculated to machine precision 
by exploiting the refinement relation of scaling functions and then tabulated for future use. 
The necessary convolution to obtain the Hartree potential requires FFTs on a grid that is 
doubled in each dimension to avoid spurious periodic interactions, but this can be performed 
without additional computational effort by modifying the FFT algorithm to exploit the fact 
that the charge density is zero on the additional grid points. This latter optimisation would 
also benefit the Cutoff Coulomb approach. A representation of the Hartree potential arising 
from the MIC Coulomb potential results that is essentially exact. 
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V. OPEN BOUNDARY CONDITIONS 



The final possibility we will consider is to change not the form of the interactions, but that 
of the boundary conditions. A careful recasting of the electrostatic terms in the Kohn-Sham 
energy functional allows us to use a form suitable for calculation with Open Boundary 
Conditions (OBCs). This is achieved by replacing the reciprocal-space evaluation of the 
core-core, Hartree and local pseudopotential energy terms by calculations performed in real 
space, which assume no periodicity of the system. 

The core-core interaction energy is calculated as a Coulombic sum of the interaction ener- 
gies of point charges as in Eq. ([s]). We describe in Appendix B how the local pseudopotential 
Viocps (i") can be calculated in real space. 

The Hartree potential VH(r) is obtained by solving the Poisson Eq. ^ in real space. 
The multigrid methocP^ represents an efficient approach for solving for the potential, given 
the charge density sampled on a regular grid and Dirichlet boundary conditions on the 
faces on the simulation cell, dQ. By using a hierarchy of successively coarser grids along 
with interpolation and restriction operators to transfer the problem between the grids, the 
multigrid approach addresses the problem of critical slowing down that plagues stationary 
iterative methodd^. For a more thorough discussion of the approach the reader is referred 
to Refs.'^^'^. In the simplest approach, second-order finite differences (FDs) are used to 
approximate the Laplacian in Eq. However, there is evidencd^^'^ that this is not suffi- 
ciently accurate for DFT calculations. One way to assess the accuracy of the solution is by 
comparing the values of two expressions for the Hartree energy, namely 



V}i{r)n{r)dr 



and 



in 
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The relative discretization error, defined as 

d 



(26) 



(27) 



(28) 



can then serve as a measure of the inaccuracy of the solution. Figure |3] shows how this error 
is unacceptably large when a second-order solver is used. The problem can be addressed 
by employing high-order defect correction, where higher-order finite differences are used 



15 



to iteratively correct the solution obtained with a second-order solvei'^. In this way the 
discretization error can be systematically reduced (Figure |3]) with moderate computational 
cost. No changes to the second-order solver are necessary. The computational cost of the 
multigrid approach scales linearly with the volume of the simulation cell, albeit with a large 
prefactor. 




Figure 3. Relative discretization error Eq. (28) in the Hartree energy vs. the order of the finite 



differences used in the defect correction of the second-order solution, on the example of aspartate. 
An order of 2 corresponds to the uncorrected solution. Smeared ions were used. 



The multigrid method does not rely on any particular form of the Dirichlet boundary 
conditions specified on dQ, however, to obtain a potential consistent with the OBCs used 
for the remaining energy terms, these should be 



^H(r) 



rdr' for r G dfl 



(29) 



Although the evaluation of the boundary conditions is straightforward, it is computationally 
costly, scaling unfavourably as 0{L'^V), which, for localised charge, implies 0{L'^N). To 
ameliorate this problem, a suitable coarse-grained approximation can be used instead of 
n{r'). Combined with evaluating Eq. (29) only for a subset of points in dQ and using 



interpolation in between, this leads to a reduction of the computational effort by 3-4 orders 
of magnitude, which brings it into the realm of feasibility. 
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The smeared-ion formalism^', where the total energy is rewritten by adding and sub- 
sequently subtracting Gaussian charge distributions centred on the cores, can be used in 
conjunction with the multigrid approach. In this case, the Poisson equation, Eq. ([6]), is 
solved for the electrostatic potential generated by the total charge density (due to electrons 
and smeared ions). As the cores neutralize a significant fraction of the electronic charge, 
the magnitude of the relevant quantities (charge density, potential) is smaller. Assuming 
the relative error incurred by the multigrid remains the same, this has the advantage of 
reducing the absolute error. The use of smeared ions, however, introduces approximations 
of its own. For a more detailed discussion of smeared ions the reader is referred to Ref.l^. 
We shall evaluate the approach with and without smeared ions. 

VI. CONVERGENCE PROPERTIES 
A. Small Molecular Systems 

We test these methods first on small-scale, simple systems to demonstrate their equiva- 
lence in the limit that all relevant parameters are accurately converged. For this, we select 
a test set of small ions molecules: a phosphate ion (P04'^~), pyridinium (CsNHg)^^, the 
amino-acid salt aspartate with a charge of — le, and the amino acid lysine with a charge 
of +le, the neutral molecules water (H2O) and potassium chloride (KCl). In this set, we 
have thus included two cations, two anions and neutral molecules with a relatively low and 
a very high dipole moment, respectively. Clearly, these small molecules are unchallenging 
calculations for linear-scaling DFT, of a size below the onset of any linear-scaling behaviour, 
but they serve to illustrate the main convergence issues in a controllable way, since it is 
here possible to make the simulation cell very much larger than the molecule, within feasible 
computational memory requirements. Illustrations are shown in Figure |4] of this test set. 

Each molecule was simulated in a cubic simulation cell initially of size 32.5ao, with a 
grid spacing of 0.5ao (equivalent to an energy cutoff of around 827eV), and with all NGWF 
radii set to J.OaQ. The density kernel was not truncated (all elements allowed to be nonzero) 
as the systems are too small for meaningful truncation. Norm-conserving pseudopotentials 
are employed for all the ions included here, and exchange-correlation is described by the 
PBE functional. We choose, throughout this work, to examine the convergence of the total 
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Figure 4. (Color online) Small molecules for initial tests, covering anions and cations and species 
with dipole moments. From left to right: phosphate, pyridinium, aspartate, lysine, potassium 
chloride and water. 
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Figure 5. Convergence of total energy with simulation cell size for a monopolar system (PO^^, 
left) and a dipolar system (KCl, right), showing the uncorrected results (black), and different 
forms of Makov-Payne correction: a) red: E{L) = Eq + q'^a/2L + B/L^\ b) green: E{L) = 
Eo + A/L + B/L^ + C/L^] c) blue: E{L) = Eq + B/L^- d) orange: E{L) = Eq + B/L^ + C/L^. 

energy, because although in practice one is most often interested in a quantity derived from 
it, such as formation or binding energies, the finite size errors made in total energies due to 
monopole or higher charges cannot be expected to cancel between (for example) reactant 
and product states, so convergence of the total energy must be obtained individually for 
each system. 



B. Convergence 

First, we examine the option of extrapolation to infinite size from calculations performed 
under PBCs. The black diamonds in Figure [5] show the uncorrected total energy of each 
of the six molecular species, calculated under PBCs. According to Makov and Paynd^^, the 
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total energy as a function of box size can be expected to behave approximately as 

in a cubic simulation cell of side L, where q is the total charge, Q is the quadrupole moment, 
and a is the Madelung constant, where for cubic cells a ~ 2.837. They suggest an approxi- 
mate correction scheme based on removing the leading order L-dependent term. However, 
there are two options for going about this in practice. 

Direct calculation of quadrupole moment Q from the density is problematic and a more 
reliable approach is to set the monopole charge q according to the known charge and then 
fit Eq and Q to data using a least-squares fit to data at multiple values of L. This produces 
the dotted black line in each figure. Alternatively, one could take into account that for a cell 
containing a molecule which is to some extent extended and may be somewhat polarisable, 
the mean dielectric constant is not equal to precisely unity. One could therefore also allow the 
coefficient of the 1/L term to vary freely, and allow an 0{L^^) coefficient as well. Examining 
Figure [5] we see that the finite size error for those species with a monopole charge follows 
E{L) = Eq + 0{L~^) fairly well as expected. The species with only a dipole moment (not a 
monopole) display a much weaker effect, which behaves as E{L) = Eq + 0{L~^). However, as 
the charge distribution varies with L, and so the coefficient Q in Eq. ( [30| ) depends weakly on 
L, the fit to the Makov-Payne form is not exact. Nevertheless, in the small charged molecules 
used here, the fitted Makov-Payne correction achieves a fairly well-defined correction to the 
total energy, aligning each individual energy to the extrapolated infinite-cell-size limit even 
for smaller cells, producing a good fit. However, the extra freedom allowed by varying q or 
introducing 0{L~^) terms are seen to produce a less useful extrapolation, by fitting to noise. 
This can only be seen for sure by comparing to the known answer obtained under one of the 
correction schemes as seen below. 

The effect of self-consistency in these small systems is not very strong: that is to say, the 
rearrangement of the charge due to the influence of the potential from neighbouring images 
of the cell is not very great. Henceforth, for Makov-Payne results, we will show the corrected 
result i?PBc(-^) ~ EyiY-iL) + Eq, where Emf{L) is the appropriate Makov-Payne choice, as 
this result falls on a comparable scale to the results for the other schemes, enabling visual 
comparison. 

Within the cutoff Coulomb approach, we can individually vary the size of the original 
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cell, the size of the padded cell, and the cutoff radius of the interaction. We note that the 
results obtained are converged to less than 1/ieV/atom once the extent of the density of 
the molecule is less than that of the original cell. Since the interaction is constructed in 
reciprocal space but has a sharp cutoff in real space at i?c, care must be taken to include 
sufficient padding that the cutoff falls within a vacuum region, and the region of 'ringing' 
induced by the cutoff is at least 5-lOao from any significant values of nonzero density. Once 
this is achieved, residual variation of the result with Rc is well below lyueV/atom. 

For the MIC approach, we obtain near-identical results for our implementation of the 
Marty na-Tuckerman approach as compared to our implementation of the approach of Gen- 
ovese et. al. We thus only show the MT approach henceforth, which was rather easier to 
parallelise in the current methodology, even though the Genovese et. al. approach is tech- 
nically more sophisticated and requires less computational effort overall due to the lower 
padding requirements. Martyna and Tuckerman note that to obtain accurate reciprocal- 
space representation of the MIC Coulomb potential, a smaller grid spacing is sometimes 
required compared to the requirements of a comparable PBC calculation. Alternatively, one 
can represent just the density and potential on a finer grid. Taking the latter approach, we 
compared grid spacings 2. Ox, 2.5x and 3. Ox the underlying psinc grid for representation 
of the density and potential. While the results do show minor variations (from 20 to 100 
/xeV/atom depending on the system), this variation is present to the same extent also in PBC 
calculations so should not be attributed to the MT approach itself — rather it is thought 
to result from changing the grid in discrete evaluation of the XC energy integral. We thus 
employ the standard 2.0 x fine grid spacing throughout the rest of this work. 

We show in Figure [6]the total energy of the test systems evaluated all the above methods. 



The CC results use the spherical cutoff of Eq. (11 ). The results for the CC and MIC methods 



converge rapidly with system size to effectively the same value. In very small simulation 
cells, below 42ao, the extent of the 'FFT box' — and thus the total extent of the charge 
distribution — is the same as that of the simulation cell. In such cases, the simulation 
cell contains very small contributions to the total electron density that wrap through the 
periodic boundaries. Therefore, even the correction schemes do not fully account for the 
absense of periodic interactions, and a quite strong dependence on L at very small L is 
seen. However, as soon as the simulation cell is large enough that the density is contained 
fully within one cell, the result is beyond that point entirely converged with system size and 
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Figure 6. Convergence of total energy with simulation cell size for test molecules, using: Cutoff 
Coulomb (green circles), Martyna-Tuckerman (orange triangles), OBCs (blue diamonds, filled when 
smeared ions were used), and MP-corrected PBCs (red squares). CC and MT results rapidly 
approach the same converged answer once the size of the cell is greater than the extent of the 
density. This converged result agrees well with the trend of the MP-corrected lines. The OBC 
results are offset by a constant due to the approximations involved in the smeared-ion representation 
and by an error whose magnitude increases with the box size (particularly evident for KCI) as a 
consequence of the numerical evaluation of the real- space pseudopotential (see Appendix B). 

independent of L. 

However, the OBC calculations evidently produce results of a somewhat lower accuracy. 
For these results, several distinct sources of inaccuracy can be distinguished. First and 
foremost, the calculation of the local pseudopotential under OBCs is performed numerically 
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and the associated error increases with the size of the simulation celL The reasons behind 
this are explained in detail in Appendix B. For the systems and box sizes encountered here, 
the magnitude of this error is 20 — 200/ieV/atom, thus it is only apparent in the plots for 
KCl, where the magnification is the highest. Second, the use of a multigrid approach to solve 
Eq. ^ introduces a discretization error. The magnitude of this error, however, can be easily 
made negligible by employing high-order defect correction, and introducing smeared ions, 
as explained earlier in section |V} cf. Fig. [3j Third, there are approximations involved in the 



generation of boundary conditions Eq. (29) for the solution of Eq (|6j). In our implementation 
we coarse-grain charge densities (electronic when smeared ions are not used, or total when 
using smeared ions) represented on a grid by replacing cubic blocks of p x p x p gridpoints 
with a single point charge located at the centre of charge of the block (thus, in general, not 



at a gridpoint). This is only done when evaluating the integral in Eq (29). With p = 5 (used 
throughout this work) the the prefactor for the calculation of the boundary conditions is 
reduced 125-fold, whereas the associated error in the energy was less than 75/ieV/atom in 
the worst case (P04^~ in the smallest box) and diminished quickly with increasing box sizes. 
For neutral systems, even with high dipole moments, this error was less than 6/xeV/atom, 
again quickly diminishing with the box size. 

Finally, the introduction of smeared ions^'^ also affects the obtained energies, as evidenced 
by the near-constant shifts between the results with and without smeared ions, observed in 
the plots. The error incurred by using smeared ions is due to the fact that certain terms in the 
formalism (e.g. the self-interaction of every smeared ion) are calculated analytically, whereas 
other terms (e.g. the local pseudopotential energy) are calculated by integrating the relevant 
quantities on a grid. Thus, the terms that are meant to cancel only do so in the limit of an 
infinitely fine grid. For the systems discussed here, the residual error is 100 — 300/ieV/atom, 
outweighing the reduction in the other sources of error that smeared ions bring about - it 
is apparent from the plots that the calculations would be more accurate without smeared 
ions. Smeared ions find use in the context of implicit solvent calculation^, as they allow 
the dielectric continuum to polarise in response not only to the electronic density, but also 
to the density of the smeared cores. For calculations in vacuum involving the systems of 
interest here, their introduction negatively impacted accuracy. 

Overall, one can conclude from these tests that in small systems, both the CC and MIC 
methods can be used with confidence, once the system size is large enough that the charge 
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density is fully contained within the appropriate box. Extrapolation-based techniques can 
correct energies to comparable accuracy, but should be used with care and the use of excessive 
variational freedom in the parameters tends to worsen results. Finally, when using OBCs, 
the energy is actually expected to very slowly diverge with the size of the simulation cell, 
due to the inaccuracies involved in the evaluation of the local pseudopotential. This effect, 
compounded by the near-constant shift due to the use of smeared ions means that OBC 
results should only be compared against other OBC results rather than results from PBC 
calculations. 

VII. LARGE SYSTEMS AND COMPUTATIONAL OVERHEAD 

To validate and compare these methods in a more realistic setting, it is necessary to 
examine their performance in larger-scale systems more typical of the real applications of 
linear-scaling DFT. These will often behave very differently from very small systems, because 
it is usually impossible to perform the calculations in a simulation cell where the dimensions 
of the cell greatly exceed that of the molecule or nanostructure. Furthermore, the scaling of 
the computational effort with system size may be very different as the balance of time spent 
in different parts of the calculation changes with the number of atoms. 

To demonstrate the accuracy of the methods, and the computational overhead and the 
scaling of each of these approaches, we have simulated two fairly large systems, each com- 
prising around 1250 atoms, which for these systems is well above the threshold at which 
linear-scaling methods offer a clear advantage in terms of reduced computational effort over 
comparable traditional DFT approaches. These systems are: a fragment of the L99A/M102Q 
mutant of the T4 Lysozyme proteirP^'^, and a nanocrystal of gallium arsenide in the wurtzite 
structure, with hydrogen termination!^. Figure [t] illustrates these systems. 

The protein fragment has a high net charge of +7e as a result of the protonation state of 
its residues at normal pH, and hence displays a strong finite size effect on the total energy if 
periodic boundary conditions are employed. This makes it difficult to calculate meaningful 
binding energies of ligands to its polar binding site. The distribution of the net charge is 
largely determined by the functional groups included and to a great extent it can be viewed 
as localised on these groups, so it is not expected to depend strongly on the system size: to 
a reasonable approximation we can treat the density of this system as fixed when we vary 
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Figure 7. Illustration of test set of large systems (a) 1234-atom fragment of the L99A/M102Q 
mutant of the T4 Lysozyme protein (b) Wurtzite-structure GaAs Nanorod of 1284 atoms, with 
hydrogen atoms terminating dangling bonds. 

the size of the simulation cell. 

The GaAs nanorod, on the other hand, has no net charge, but the underyling wurtzite 
structure, with no inversion symmetry, means that when truncated at each end of the rod 
along the c-axis, Ga and As faces are exposed on opposite ends of the rod. No matter how 
the surface is terminated (in the case studied here, all dangling bonds are saturated with 
hydrogen), there will be some form of charge transfer between the ends and a dipole moment 
along the c-axis will result. If such a rod is simulated in a box of size comparable to the rod 
itself under PBCs, then the rod is effectively surrounded by an inifinite array of replicas, 
producing a very different electric field from that of an isolated rod. Indeed, unless the box 
is very large along all axes, the Ga-terminated end of the rod will be in closer proximity to 
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the As-terminated ends of rods in adjacent cells than to the As-terminated end on the on 
the same rod (and vice versa), and the rod is strongly polarisable. This is clearly a very 
different situation from the ideal situation many correction methods assume, of a strongly 
localised fixed charge distribution in a box considerably larger than the charge distribution 
itself. Because the magnitude of the dipole moment depends sensitively on the balance of 
charge distribution and the density of states at the polar surfaces of the rod, its value can 
be affected by the field created by the charge distribution of periodic images of the rod, 
bringing self-consistent effects into play. 

To perform these large simulations, we use in both cases a grid spacing of O.Sao; equiv- 
alent to a plane-wave cutoff of around 827eV, and standard well-tested norm-conserving 
pseudopotentials for each species. For the protein system, four NGWFs of radius 7.0ao were 
placed on each C, N, O and S atom, and one on each H atom. For the nanorod, larger NG- 
WFs were required to achieve good convergence, so = lOao was used, with four NGWFs 
on Ga and As and one on H. 

The total energy of the the protein fragment as a function of supercell side length is 
shown in Figure |8| We use a series of cells up to L = 400ao in size, so as to be able to 
accurately extrapolate to L — )■ oo. We see that on the larger scale (top figure), the results 
for all three methods lie on apparently the same line, agreeing with the extrapolation of the 
Makov-Payne form to large L. However, zooming in reveals two significant details: firstly, 
there remains considerable residual variation in the Makov-Payne corrected results, which 
do not converge to better than O.OSeV until L = 200ao. When they do so, they agree well 
with the MP extrapolation. The OBC result suffers from considerable residual error, mostly 
due to the approximations involved in the evaluation of the local pseudopotential, which for 
the smaller box sizes cancels out, to a degree, with the error due to the approximations in 
the construction of the boundary conditions, but for larger boxes causes the energy to very 
slowly diverge. The almost constant shift in energy incurred by the use of smeared ions is 
approximately 400yueV/atom. The CC and MT results agree very well with each other, to 
around the ImeV level. We conclude that for monopolar systems with an approximately 
fixed charge distribution, the CC and MT methods are both reliable as long as the cell is 
made large enough for the conditions of each method to hold. 

The total energy of the the nanorod as a function of supercell side length is shown in 
Figure [9j Here the default supercell is not chosen to be cubic as this would be particularly 
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Figure 8. Convergence of total energy with cell size for T4 lysozyme fragment, showing results for 
Cutoff Coulomb (green circles), Minimum Image Convention (orange triangles). Open Boundary 
Conditions (blue diamonds, filled when smeared ions were used) and Periodic Boundary Conditions 
(red squares) corrected using the Makov-Payne form (a) fitted by least-squares fitting to £"0 and B. 

inefficient for such a high-aspect ratio nanostructure. We start with = 240ao, Ly, = 
65ao as the smaUest cell able to completely enclose the rod with around lOoo padding in 
all directions, and then Ly and Lz are scaled commensurately with L^. We see that in this 
case, with a highly polarisable rod, the same Makov-Payne fit that successfully described 



the dipolar systems in Section VI now fails significantly for all the cells studied here and 
returns an Eq which does not match the L — t- 00 limit, nor does it match the CC or MIC 
results. The latter are well-converged with respect to L^, and are in good agreement with 
each other. However, again the OBC results are strongly size-dependent, as a result of the 
approximations made in order to obtain feasible computational effort at this large scale. 
The validity of the convergence of the approximate methods starts to break down beyond 
Lx ~ 300ao, resulting in significant error. 

By examining the behaviour of the dipole moment of the rod along its length, calculated 
as dx = J^xn{r)dr, we see immediately why this is the case: the dipole moment varies 



26 



-200612.5 
-200613.0 



> -200613.5 



-200614.0 
-200614.5 



II 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

O ■ A 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [ 1 1 1 1 1 1 III 

■ ■ ■ - ' 

A S J 


' ■ 




♦ 

■ o 

: oo 




iilllillllllllll|i|illlllllillllB|& 


' o 

' ■ 1 ■ 1 1 t 1 t 1 1 1 . 1 1 t ■ 1 ■ 1 1 


GaAs nanorod : 



200 300 400 500 600 700 800 900 1000 
Cell side, (a.u.) 



Figure 9. Convergence of total energy with cell size for GaAs Nanorod, showing results for Cutoff 
Coulomb (green circles), Minimum Image Convention (orange triangles). Open Boundary Condi- 
tions (blue diamonds, filled when smeared ions were used) and Periodic Boundary Conditions (red 
squares) corrected using the Makov-Payne form (c) fitted by least-squares fitting to Eq and B. 
CC results are independent of cell size for all Lx greater than the rod length. The MIC approach 
requires > 2 x -Lrodj so results are only shown for > 480ao. 

strongly with cell size because of the induced polarisation caused by the periodic images, 



as seen in Figure 10 The periodic images of the nanorods are all aligned, so if the rods 
are very close end-to-end they will tend to increase the dipole moment. However, if they 
are closer side-to-side the dipole field of the periodic image (in the opposite direction to the 
polarisation, as viewed outside the rod to its side) will tend to depolarise the rod and the 
dipole moment will decrease. Therefore there is a strong dependence of on both and 
Ly, Lz- Both of these are spurious effects when one wishes to simulate an isolated rod. We 
see that all three approaches, CC, MIC, and OBC, all correct these influences and obtain 
the 'correct' isolated result for dx even for small system sizes. 

We therefore conclude that in such cases of large, polarisable systems with a strong dipole 
moment, there is no choice but to use an approach including the truncation of periodic im- 
ages: in analogy to the study of polar thin film^, a correction scheme must be employed if 
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Figure 10. Dipole moment dx (see text) as a function of cell size Lx for a GaAs nanorod. Here 
the inset illustration is shown approximately to-scale with the L-axis. The exact form of dx{L) 
would depend on aspect ratio, and would be difficult to accurately extrapolate to L — )• oo. The 
cutoff Coulomb and MIC approaches obtain converged results for all cell sizes large enough for their 
methods to hold, while the periodic results converge only very slowly to this isolated value. 



reliable results are to be obtained. We have demonstrated that the approaches of Coulomb 
truncation and MIC are suitable for this purpose. The inaccuracies inherent in the OBC 
approach are particularly visible in the case of the nanorod, as the very large box sizes cause 
the error associated with the evaluation of the local pseudopotential to become unaccept- 
ably large. The error due to the smeared ion representation is approximately 700//eV/atom 
and for the smallest box sizes it conveniently cancels most of the error in the local pseu- 
dopotential, however for the larger simulation cells the energy inevitably diverges. Although 
this divergence is slow (compared to the total energy of the system), in the absence of a 
monopole charge and the associated 0(L~^) term it makes the OBC results unacceptably 



inaccurate for energies. Figure 10 shows nevertheless that for other physical properties, such 



as the dipole moment, it may be reliable. 
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VIII. CONCLUSIONS 



We have described and applied three different methods, each with a rather different 
theoretical basis, to the study of calculations on charged and dipolar systems using hnear- 
scaling density functional theory under periodic boundary conditions. We have shown that 
with each of these methods it is possible, while staying within a nominally periodic formalism, 
to achieve the desired limit of equivalence of any calculated properties to those of a single 
isolated system. 

In small systems, post-hoc correction schemes are capable of extrapolating to the isolated 
hmit on the basis of several calculations performed under PBCs, but only if simulation cells 
are used which arc very large compared to the system being studied. The first-order term 
of the Makov-Payne correction, on its own, is inadequate for accurate results, but by fitting 
the coefficient of the 0{L~^) term, one can acheive an accurate result for a cubic cell as long 
as there is not a dipole moment present of comparable physical size to the cell itself. This 
is clearly a highly computational expensive approach due to the need to simulate several 
large cells to achieve an accurate fit, and is not really practical for production calculations. 
Fitting further coefficients of the Makov-Payne expansion tends to reduce the accuracy by 
over-fitting to numerical noise. 

However, we have also seen that the larger systems encountered in large linear-scaling 
DFT calculations can behave very differently from the small molecules in the test set. In 
particular, there is scope in large systems for considerable long-range charge redistribution 

in response to the effect of periodic images, so reliable extrapolation from simulations using 
a small unit cell are then impossible. In such situations, one has no choice but to use a 
method that explicitly negates the effect of periodic images. 

Approaches which redefine the Coulomb potential to avoid periodic interactions, either 
by using the Minimum Image Convention (whether in the form proposed by Martyna and 
Tuckerman, or in the rather different form by Genovesc et al.) or the cutoff Coulomb 
method, have been seen to rapidly converge to the isolated result as soon as the conditions 
required as outlined above are met. In the case of the MT formulation, this is that the size 
of the simulation cell be at least twice the extent of the electron density in a given direction, 
while in the Genovese form, this requirement is relaxed due to the method being performed 
on what is effectively a padded real-space grid. 

29 



The cutoff Coulomb approacfi is seen to produce accurately converged results for a single- 
shot calculation, regardless of the size of the simulation cell (as long as it is bigger than the 
extent of the nonzero density). The only requirement is that the original cell must, for 
the purposes of Fourier transforms, be embedded in a padded cell of sufficient size. This 
generally entails quite a large temporary memory requirement, and in small systems the 
performance of FFTs can become the limiting factor on the speed of the whole calculation. 
However, in large systems, where the Hartree calculation is generally not a significant part 
of the total computational effort, this is no longer an issue. 

Finally, the use of Open Boundary Conditions, while exact in principle, is seen to entail 
several further approximations in practice to render it feasible. In particular, these enter 
into the evaluation of the Dirichlet boundary values on the faces of the simulation cell, and 
the use of a smeared-ion representation and the evaluation of the local pseudopotential in 
real space. These approximations combine to give a residual finite-size error notably larger 
than the other methods can achieve. Furthermore, the multigrid approach to the Hartree 
potential is computationally rather demanding and does not parallelise as well as the rest 
of the approach. This makes the OBC approach the least efficient method presented here. 
However, it has one major advantage the others cannot match, namely that it can be used 
with an nonhomogeneous dielectric constant, in the context of implicit solvent calculations. 
For future calculations of this type, further investigation will be required in order to develop 
means to calculate the boundary conditions to higher accuracy - such as by combining the 
multigrid OBC approach with one of the other schemes here solely for the determination of 
boundary conditions. 

We noted also that two of the methods considered here can benefit from similar speedups 
by suitable treatment of Fourier transforms padded with zeroes. In both the cutoff Coulomb 
approach and the MIC approach, there is a need to perform a FFT of the charge density to 
reciprocal space under conditions where the value on most of the real-space grid is known 
to be zero. In such cases, it has been shown that algorithms can be designed which not 
only significantly reduce the computational expense of such a transform but also reduce 
the memory usage by not explicitly storing the zeros. The MIC implementation of Gen- 
ovese et al. employs such an approach, but the Martyna-Tuckerman and cutoff Coulomb 
approaches could in principle also do so. This would render them all very similar in terms 
of computational cost, only marginally above that of the original, uncorrected approach. 
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APPENDIX A: FOURIER COEFFICIENTS OF THE 
CYLINDRICALLY-CUTOFF INTERACTION 

The integral for the Fourier components vcc{G) of the Coulomb interaction cut off over 
a cylinder of length 2L and radius R can be written 




d(j)dxdp . 



Here we have taken the cylinder to be aligned along x, and taken Gp to lie in the xy-plane, 
without loss of generality. To ensure that the resulting expression is finite and well-behaved 
for all non-negative values of G, we identify four regions which must be treated separately: 

Gp >0,G^> 0; 
Gp = 0, > 0; 
Gp > 0, = 0; 
Gp = 0, = 0. 

The latter three cases all allow significant simplification of the integral and will be examined 
first. 

The Gp — 0, Gx — terms are the only ones where the integral can be performed fully 
analytically: 
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The Gp = 0, Gx > terms can be rendered into a well-behaved integral over x: 



vcc{G) 



R rL 



2^ p^iG^x 



10 J-LJo ^/p'^ + x^ 

-L 



d0da;dp 



_^ 



= 47r^ (yR 



^ + x'^ — x^ COS GxX dx 



which can be evaluated numerically with no significant difiiculties 

Similarly, the Gp > 0, G^ — terms can be made into a well-behaved integral over p: 
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R rL r'i 
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pJo{Gpp) dp. 



which also remains well-behaved over its range. 

Finally, for Gp > 0, > 0, the integral cannot so easily be put in a 1-dimensional form 
for easy evaluation. However, if the cylinder length L is first taken to infinity (effectively 
making the interaction periodic in x), the integrals become tractable, then the resulting 
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answer can be convolved with a top-hat function to retrieve the desired hmits on the integral. 
The top hat function is defined in terms of the Heaviside step function: 

T(r) = Q{x + L) - Q{x - L) . 

The transform of the Coulomb interaction for the infinite cylinder would give 

nR r-oo /•2n iGpp sin (t>+iGxX 

vic{G)= / / ^^—==^ d<f>dxdp , 



so we can write the transform of the finite cylinder as 

f'H noo n2n 

vcc{G)^ / / / 2'(r)^;/c(r)e*(^'"'''°*+^^^M0dxpdp. 

Jo J-oo Jo 

By the convolution theorem we can write the transform of the product of two functions in 
real space as the convolution of these two functions in reciprocal space. Using H for our 
primed set of reciprocal space coordinates we get: 

^^^^^^ ^ (2^ / ^^^(^)^(^ - ^) ■ 
All three integrals for vic{ii) can be done analytically: 

rR rco p2tv 

v:c{ii) = III - cos{Hxx) cos{Hpp sin 0) d^dxdp 

Jo J-odJo Vp^ + x^ 



2 / pKo{H^p)cos{Hppsm(f))d(f)dp 
Jo Jo 

= 47r /" p Ko{H^p) Jo (Hpp) dp 
Jo 



= 47r 



1 + HpRKo{H,R) Ji{HpR) - H,RK,{H,R) UH^R) 



Hj + HI 



This expression is in fact very simple to evaluate as it contains no Bessel functions of 
higher order than 1. These can be rapidly evaluated using accurate polynomial approxima- 
tions over the domain required for the integrals. 

For the step function, the transform is well known 

f{G)^j eicp[iG^x]dxS{Gp) 
2sin(G^L) 



Combining the two gives us 



-S{G,) . 



After performing the Hp integral to leave only Hp — Gp, we obtain 
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Only the latter integral term needs to be calculated numerically. One can see that as 
R ^ oo and L — )■ oo, the modified Bessel function terms tend to zero, leaving only the 
expected Att/G"^ behaviour from the first part. When performing the integral numerically, 
the denominator damps out the oscillations rapidly so the region of integration can be 
relatively small. A fairly fine mesh must be used to capture the oscillations of the sine 
function, but not unmanageably so for the G-vectors typically required. We used 200001 
points in this work, ensuring convergence to 10 significant figures for the largest Gx and L 
values required. 



APPENDIX B: CALCULATION OF THE LOCAL PSEUDOPOTENTIAL IN 
REAL SPACE 

The local pseudopotential VJocps (r) can be evaluated in real space as a sum of spherically- 
symmetrical contributions from all atomic cores /, each located at Rf. 

Mocps (r) = ^-P«.^ (1^ - ^^1) ■ (31) 

To generate the local pseudopotential Viocps,/ (r) due to core / at a point r in real space, the 
continuous Fourier transform can be employed: 

Sloops,/ (r - Ri) = [ Viocps,/ (G) e'^ ('-^^)dG 

(27r) J 

[Zn) J 
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where we have set x = r — R/. We then use the expansion of the plane wave e'*^ '' in terms 
of locahsed functions, to obtain: 



Vlocps,/ (x) - 3 / Viocps,/ (G) 

(2vr) J 



I 



dG 



X 47r ^ ^ i^ji (Gx) Zim (^g) Zim i^^ 

1=0 m=~l 

where ji are spherical Bessel functions of the first kind and Zi^ are the real spherical har- 
monics. A simple rearrangement leads to 



Vfocps,/ (x) = 3 N N i Zlm {^y 

(27r) i=Om=-l 



X J V^ocps,/ (G) ji (Gx) Zlm {^g) dG . (33) 

After changing into spherical polar coordinates and applying the orthonormality property 
of spherical harmonics, the above expression simplifies to a spherically-symmetric form: 

oo 

^ocps,/ (x) = [ V^ocps,/ (G) dG . (34) 

(ztt) J X 







In practice, it is sufficient to evaluate this expression once, for every ionic species s{I), 
rather than for every core /, on a fine radial grid with x ranging from to a maximum value 
dictated by the size of the simulation cell in use. A finite upper limit, Gcut, corresponding 
to the longest vector representable on the reciprocal grid, should be used in the integral in 



Eq (34), in order to avoid aliasing when transforming from reciprocal to real space. 



The numerical evaluation of the integral in Eq. (34) is not straightforward. One source 
of difficulties is the oscillatory nature of sin(Gx). For larger cells, the oscillations become 
so rapid that the resolution with which the reciprocal-space coefficients Viocps.s (G) of the 
pseudopotential are provided, typically 0.05 A~^, is not sufficient and it becomes necessary to 
interpolate Viocps,s (G), and the whole integrand, in between these points. Another difficulty 
is caused by the singularity in Viocps.s {G) as G — t- 0, where the behaviour of Viocps.s (G) 
approaches that of —Z^/G"^ (where Zg is the charge of the core of species s). Although 
the integral is convergent, this singularity cannot be numerically integrated in an accurate 
fashion, and it also contributes to making the abovementioned interpolation inaccurate 
at low G's. This is partially alleviated by subtracting the Coulombic potential, —Zg/G"^, 
before interpolating to the fine radial reciprocal-space grid, and then adding it back but 
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the residual numerical inaccuracy leads to a near-constant shift of the obtained real-space 
pseudopotential, which in turn results in errors in the total energy in the order 0.01%. 
To address this problem the pseudopotential can be partitioned into a short-range and a 



long-range term, similarly as in Eq. (15). This leads to 



oo 

( \ _ erf {ax) An f ~ 
nocps,/ \^) — 1- , ,3 / t/iocps,/ Kyj) 







1 - exp ( — ^ 



^^"(^"^-GrfG, (35) 



X 



4a2 

where the first term, with the error function, is the long range part and the second is the 
short range part and a is an adjustable parameter that controls where the transition between 
short-range and long-range takes place. 



Owing to the 



1 - exp 



factor, the singularity at G = is avoided in the same 
way as in Eq. (20) and the integral can be accurately evaluated numerically, provided a 
is large enough. Smaller values of a make the numerical integration less accurate, because 
the oscillations at low values of G increase in magnitude. Larger values of a increase the 
accuracy of the integration, however, they lead to a faster decay of the reciprocal-space term 
and cause the long-range behaviour to be increasingly more dictated by the first term in 



the RHS of Eq. (35). As this term is calculated in real space, it lacks the oscillations that 
are expected to be present in the pseudpotential at large x, due to the finite value of Gcut, 
causing aliasing. For this reason a needs to be as small as possible, without negatively 
impacting on the accuracy of the numerical integration. 

The accuracy of the approach can be assessed by comparing the real-space tail of the 
obtained pseudopotential with the Coulombic potential. Since the obtained pseudopoten- 
tial is expected to oscillate slightly so that it takes values above and below —Zg/x, a good 

measure of accuracy, which we will call 6, is ^ locps.s ^' -^—^1 where the average 

runs over the real-space tail of the pseudopotential, from, say, 5ao to the maximum x for 
which V^ocps.s {x) is evaluated. Ideally, h should be zero. Numerical inaccuracy will cause a 
shift in Viocps,s {x) which will present itself as a finite, non-zero value of h. Direct numerical 



integration of Eq. (34) using various high order quadrature schemes results in values of h in 
the order of 0.01, which can be reduced by an order of magnitude by interpolating to a very 
fine radial reciprocal-space grid. Subtracting the Coulombic potential and integrating only 
the difference between Viocps,s (G) and the Coulombic potential numerically, while analyti- 
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cally integrating the remaining part reduces b to about 0.0005. Application of the proposed 
approach Eq. (35) yields 6 = 5 x 10^^ for a = 0.5// and b = 3 x 10^^ for a = 0.1//, where 
/ is the length of the simulation cell. The total energy is then insensitive (to more than 
0.0001%) to the choice of a, provided it is in a wide "reasonable" range of 0.1// — 2/1. 
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